---
title: "Figure 1"
output: 
---

# Figure 1

# Who's to Blame? Postconflict Violence and Public Attitudes Towards Peace Agreements
# Wyer, Frank. 

#clear environment
```{r clear environment}
rm(list = ls())
```

# uncomment and set working directory to replication archive
# setwd("~/blame_replication")

# Uncomment to install packages if necessary
# install.packages("tidyverse")
# install.packages("estimatr")

#load packages
```{r}
library(tidyverse)
library(estimatr)
```

#read in data
```{r}
survey_clean <- read.csv("survey_clean.csv")
```

#estimate models for treatment effects versus control condition
```{r estimate treatment effects vs control}
eln_shortmodel_95 <- lm_lin(formula = eln_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
eln_shortmodel_90 <- lm_lin(formula = eln_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

dissident_shortmodel_95 <- lm_lin(formula = dissident_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
dissident_shortmodel_90 <- lm_lin(formula = dissident_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

accords_shortmodel_95 <- lm_lin(formula = accords_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
accords_shortmodel_90 <- lm_lin(formula = accords_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

index_shortmodel_95 <- lm_lin(formula = outcomes_zscale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
index_shortmodel_90 <- lm_lin(formula = outcomes_zscale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)
```

#releveled treatment variable where postconflict violence condition (T1) is baseline instead of control
```{r treatment with T1 as baseline}
#create treatment variable with T1 as reference instead of control 
survey_clean$treatment_vsT1 <- as.factor(survey_clean$treatment)
survey_clean$treatment_vsT1 <- relevel(survey_clean$treatment_vsT1, ref = "T1")
```

#estimate models for treatment effects versus postconflict violence condition (T1)
```{r compute treatment effects vs T1}
eln_diff_95 <- lm_lin(formula = eln_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
eln_diff_90 <- lm_lin(formula = eln_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

dissidents_diff_95 <- lm_lin(formula = dissident_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)

dissidents_diff_90 <- lm_lin(formula = dissident_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

accords_diff_95 <- lm_lin(formula = accords_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)
accords_diff_90 <- lm_lin(formula = accords_scale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)

index_diff_95 <- lm_lin(formula = outcomes_zscale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)

index_diff_90 <- lm_lin(formula = outcomes_zscale ~ treatment_vsT1, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)
```

#extract treatment effects of postconflict violence treatment (T1) for dataframe
```{r extract treatment effects of postconflict violence treatment}
T1_effects <- rbind(
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_95$coefficients['treatmentT1'], Conf_Low = eln_shortmodel_95$conf.low['treatmentT1'], Conf_High = eln_shortmodel_95$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_90$coefficients['treatmentT1'], Conf_Low = eln_shortmodel_90$conf.low['treatmentT1'], Conf_High = eln_shortmodel_90$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_95$coefficients['treatmentT1'], Conf_Low = accords_shortmodel_95$conf.low['treatmentT1'], Conf_High = accords_shortmodel_95$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_90$coefficients['treatmentT1'], Conf_Low = accords_shortmodel_90$conf.low['treatmentT1'], Conf_High = accords_shortmodel_90$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_90$coefficients['treatmentT1'], Conf_Low = index_shortmodel_90$conf.low['treatmentT1'], Conf_High = index_shortmodel_90$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_95$coefficients['treatmentT1'], Conf_Low = index_shortmodel_95$conf.low['treatmentT1'], Conf_High = index_shortmodel_95$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_95$coefficients['treatmentT1'], Conf_Low = dissident_shortmodel_95$conf.low['treatmentT1'], Conf_High = dissident_shortmodel_95$conf.high['treatmentT1']), 
data.frame(Treatment = "Postconflict Violence Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_90$coefficients['treatmentT1'], Conf_Low = dissident_shortmodel_90$conf.low['treatmentT1'], Conf_High = dissident_shortmodel_90$conf.high['treatmentT1']))
```

#extract treatment effects of rebel culpability treatment (T2B) for dataframe
```{r extract treatment effects of rebel culpability treatment}
T2B_effects <- rbind(
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_95$coefficients['treatmentT2B'], Conf_Low = eln_shortmodel_95$conf.low['treatmentT2B'], Conf_High = eln_shortmodel_95$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_90$coefficients['treatmentT2B'], Conf_Low = eln_shortmodel_90$conf.low['treatmentT2B'], Conf_High = eln_shortmodel_90$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_95$coefficients['treatmentT2B'], Conf_Low = accords_shortmodel_95$conf.low['treatmentT2B'], Conf_High = accords_shortmodel_95$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_90$coefficients['treatmentT2B'], Conf_Low = accords_shortmodel_90$conf.low['treatmentT2B'], Conf_High = accords_shortmodel_90$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_90$coefficients['treatmentT2B'], Conf_Low = index_shortmodel_90$conf.low['treatmentT2B'], Conf_High = index_shortmodel_90$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_95$coefficients['treatmentT2B'], Conf_Low = index_shortmodel_95$conf.low['treatmentT2B'], Conf_High = index_shortmodel_95$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_95$coefficients['treatmentT2B'], Conf_Low = dissident_shortmodel_95$conf.low['treatmentT2B'], Conf_High = dissident_shortmodel_95$conf.high['treatmentT2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_90$coefficients['treatmentT2B'], Conf_Low = dissident_shortmodel_90$conf.low['treatmentT2B'], Conf_High = dissident_shortmodel_90$conf.high['treatmentT2B']),

data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_diff_95$coefficients['treatment_vsT1T2B'], Conf_Low = eln_diff_95$conf.low['treatment_vsT1T2B'], Conf_High = eln_diff_95$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_diff_90$coefficients['treatment_vsT1T2B'], Conf_Low = eln_diff_90$conf.low['treatment_vsT1T2B'], Conf_High = eln_diff_90$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_diff_95$coefficients['treatment_vsT1T2B'], Conf_Low = accords_diff_95$conf.low['treatment_vsT1T2B'], Conf_High = accords_diff_95$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_diff_90$coefficients['treatment_vsT1T2B'], Conf_Low = accords_diff_90$conf.low['treatment_vsT1T2B'], Conf_High = accords_diff_90$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_diff_90$coefficients['treatment_vsT1T2B'], Conf_Low = index_diff_90$conf.low['treatment_vsT1T2B'], Conf_High = index_diff_90$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_diff_95$coefficients['treatment_vsT1T2B'], Conf_Low = index_diff_95$conf.low['treatment_vsT1T2B'], Conf_High = index_diff_95$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissidents_diff_95$coefficients['treatment_vsT1T2B'], Conf_Low = dissidents_diff_95$conf.low['treatment_vsT1T2B'], Conf_High = dissidents_diff_95$conf.high['treatment_vsT1T2B']), 
data.frame(Treatment = "Rebel Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissidents_diff_90$coefficients['treatment_vsT1T2B'], Conf_Low = dissidents_diff_90$conf.low['treatment_vsT1T2B'], Conf_High = dissidents_diff_90$conf.high['treatment_vsT1T2B'])
)
```

#extract treatment effects of govt culpability treatment (T2A) for dataframe
```{r extract treatment effects of govt culpability treatment}
T2A_effects <- rbind(
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_95$coefficients['treatmentT2A'], Conf_Low = eln_shortmodel_95$conf.low['treatmentT2A'], Conf_High = eln_shortmodel_95$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_shortmodel_90$coefficients['treatmentT2A'], Conf_Low = eln_shortmodel_90$conf.low['treatmentT2A'], Conf_High = eln_shortmodel_90$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_95$coefficients['treatmentT2A'], Conf_Low = accords_shortmodel_95$conf.low['treatmentT2A'], Conf_High = accords_shortmodel_95$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_shortmodel_90$coefficients['treatmentT2A'], Conf_Low = accords_shortmodel_90$conf.low['treatmentT2A'], Conf_High = accords_shortmodel_90$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_90$coefficients['treatmentT2A'], Conf_Low = index_shortmodel_90$conf.low['treatmentT2A'], Conf_High = index_shortmodel_90$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_shortmodel_95$coefficients['treatmentT2A'], Conf_Low = index_shortmodel_95$conf.low['treatmentT2A'], Conf_High = index_shortmodel_95$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "95%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_95$coefficients['treatmentT2A'], Conf_Low = dissident_shortmodel_95$conf.low['treatmentT2A'], Conf_High = dissident_shortmodel_95$conf.high['treatmentT2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Control", Conf_Level = "90%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissident_shortmodel_90$coefficients['treatmentT2A'], Conf_Low = dissident_shortmodel_90$conf.low['treatmentT2A'], Conf_High = dissident_shortmodel_90$conf.high['treatmentT2A']),

data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_diff_95$coefficients['treatment_vsT1T2A'], Conf_Low = eln_diff_95$conf.low['treatment_vsT1T2A'], Conf_High = eln_diff_95$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Support Peace Agreement with ELN", Estimate = eln_diff_90$coefficients['treatment_vsT1T2A'], Conf_Low = eln_diff_90$conf.low['treatment_vsT1T2A'], Conf_High = eln_diff_90$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_diff_95$coefficients['treatment_vsT1T2A'], Conf_Low = accords_diff_95$conf.low['treatment_vsT1T2A'], Conf_High = accords_diff_95$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Confidence in 2016 Peace Agreement", Estimate = accords_diff_90$coefficients['treatment_vsT1T2A'], Conf_Low = accords_diff_90$conf.low['treatment_vsT1T2A'], Conf_High = accords_diff_90$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_diff_90$coefficients['treatment_vsT1T2A'], Conf_Low = index_diff_90$conf.low['treatment_vsT1T2A'], Conf_High = index_diff_90$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Peace Agreement Attitudes Index*", Estimate = index_diff_95$coefficients['treatment_vsT1T2A'], Conf_Low = index_diff_95$conf.low['treatment_vsT1T2A'], Conf_High = index_diff_95$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "95%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissidents_diff_95$coefficients['treatment_vsT1T2A'], Conf_Low = dissidents_diff_95$conf.low['treatment_vsT1T2A'], Conf_High = dissidents_diff_95$conf.high['treatment_vsT1T2A']), 
data.frame(Treatment = "Government Culpability Treatment", Versus= "vs Postconflict Violence Treatment", Conf_Level = "90%", Outcome = "Support Peace Agreement with FARC Dissidents", Estimate = dissidents_diff_90$coefficients['treatment_vsT1T2A'], Conf_Low = dissidents_diff_90$conf.low['treatment_vsT1T2A'], Conf_High = dissidents_diff_90$conf.high['treatment_vsT1T2A'])
)
```

#combine all treatment effects into single dataframe
```{r merge treatment effects into a single dataframe}
treatment_effects_all <- rbind(T1_effects, T2A_effects, T2B_effects)
```

#convert characters to factors for plot ordering
```{r make character columns factors to control order for plot}
treatment_effects_all$Outcome <- factor(treatment_effects_all$Outcome, levels = c("Peace Agreement Attitudes Index*", "Support Peace Agreement with FARC Dissidents", "Support Peace Agreement with ELN", "Confidence in 2016 Peace Agreement"))
treatment_effects_all$Treatment <- factor(treatment_effects_all$Treatment, levels = c("Postconflict Violence Treatment", "Rebel Culpability Treatment", "Government Culpability Treatment"))
treatment_effects_all$Versus <- factor(treatment_effects_all$Versus, levels = c("vs Postconflict Violence Treatment", "vs Control"))
```

#plot all treatment effects
```{r combined plot code}
ggplot(treatment_effects_all %>% pivot_wider(names_from = "Conf_Level", values_from = c(Estimate, Conf_Low, Conf_High)), aes(y = Outcome, x = `Estimate_95%`, shape = Versus)) +
        geom_vline(xintercept = 0, 
                   colour = gray(1/2), lty = 2) +
scale_shape_manual(name = "", values = c(15, 19), breaks = c("vs Control", "vs Postconflict Violence Treatment"), labels = function(x) str_wrap(x, width = 15)) +
        geom_point(aes(y = Outcome, 
                    x = `Estimate_95%`), position=position_dodge(.45), size = 2.5) + 
        geom_linerange(aes(y = Outcome, 
                     xmin = `Conf_Low_90%`,
                     xmax = `Conf_High_90%`),
                   lwd = 1, position=position_dodge(.45)) + 
        geom_linerange(aes(y = Outcome, 
                     xmin = `Conf_Low_95%`,
                     xmax = `Conf_High_95%`),
                   lwd = 1/2, position=position_dodge(.45)) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 15)) +
scale_x_continuous(breaks = c(-.6, -.4, -.2, 0, .2)) +
    guides(shape = guide_legend(byrow = TRUE)) +
                           labs(x = "", y = "", shape = "") + 
theme_bw() +
    facet_wrap(Treatment ~ ., strip.position = "top", labeller = label_wrap_gen(width = 15, multi_line = TRUE)) +
        theme(legend.spacing.y = unit(1, "cm"), strip.text = element_text(face = "bold", size = 10))
```
